Asymptotics of work distributions in non-equilibrium systems 
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The asymptotic behaviour of the work distribution in driven non-equilibrium systems is deter- 
mined using the method of optimal fluctuations. For systems described by Langevin dynamics the 
corresponding Euler-Lagrange equation together with the appropriate boundary conditions and an 
equation for the leading pre-exponential factor are derived. The method is applied to three repre- 
sentative examples and the results are used to improve the accuracy of free energy estimates based 
on the application of the Jarzynski equation. 
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I. INTRODUCTION 

Recent progress in the statistical mechanics of non- 
equihbrium systems centered around fluctuation [l|, |2| 
and work [3|, |j| theorems has profound implications for 
both theory and applications. Rather complementary 
to the traditional emphasis of statistical mechanics on 
typical behaviour of systems these new lines of research 
put the large deviation properties of thermodynamic vari- 
ables like work or entropy into focus. Of particular in- 
terest for many practical applications is the use of the 
Jarzynski equation [3| 



-/3AF 



-PW\ 



(1) 



to determine the free energy difference IS.F between two 
equilibrium states at inverse temperature [3 from the 
work distribution P{W) characterizing irreversible tran- 
sitions between these states. The method works best if 
AF is of the order of the thermal energy, 1//3. Detailed 
knowledge of free energy differences in mesoscopic sys- 
tems is of crucial importance for problems like the con- 
formations of polymers, the decay of metastable states, 
or the efficiency of molecular motors. 

It is very remarkable and particularly attractive for 
systems with long relaxation times that equilibrium in- 
formation like Ai^ may be obtained from fast changes 
of state. The method has been successfully employed in 
experiments on mesoscopic systems [a H S Q ^^ well as 
in numerical simulations 3, 10| where, however, its su- 
periority to other methods is still under debate [ll|] . The 
main problems arise from the exponential average in ([T]) 
which is dominated by small values of W from the tail of 
the distribution P[W). Since these large deviations are 
rarely sampled the resulting free energy estimate may be 
poor. An equivalent observation is |12l.ll3l| that the dom- 
inant trajectories contributing most to the average in the 
r.h.s of ID) are in general rather different from the typical 
ones, i.e. from those with the highest probability. Several 
methods have been put forward to improve the accuracy 
of free energy estimates by, e. g., i nclu ding information 
from the backward process [ J , Il4l. llSl llq . using map- 
pings and auxiliary drifts [iTI . llSl . |l9[ , or implementing 
biased path ensembles [20l l2ll . |23 ] . 



In the present paper we devise a method to analytically 
determine the asymptotics of the work distribution P{W) 
of driven Langevin systems for very small or large values 
of W . We demonstrate that fitting these asymptotics to 
the region of work values that is still sufficiently sam- 
pled by experiment or simulation significantly improved 
estimates of the free energy difference may be obtained. 



The procedure builds on the method of optimal fluctu- 
ation which rests on the general assumption of large de- 
viation theory [23, [2J| that the probability of an unlikely 
event is dominated by the most probable fluctuation giv- 
ing rise to it. All other possibilities to bring the same 
result about are even more unlikely and may be safely 
neglected. In physical context the method was originally 
proposed to determine the asymptotic tail of the elec- 
tronic density of states in random potentials [2^ |2a 123] . 
Later applications include the motion of charge density 
waves in disordered media 1281 . the velocity distribution 
in Burgers turbulence [29|, |30|, anomalous optical ab- 
sorption in disordered semiconductors [31|, and the free 
energy distribution of a directed polymer in a random 
medium [34] . Recently there have also been applications 
to optimal control theory [32] and error correcting codes 
[33] ■ In the present example of Langevin dynamics the 
method corresponds to a saddle-point approximation in 
a functional integral over stochastic trajectories. 



The paper is organized as follows. In section 2 we 
outline the general theory. First the optimal path for 
a work value in the tail of P{W) is determined by the 
solution of a variational problem, then the contribution 
from neighbouring paths is included. Section 3 concerns 
the discussion of three concrete examples. For the first 
the complete P{W) can be determined analytically so it 
serves merely as a test of our method. In the second we 
compare our results with numerical simulations of the 
Langevin equation whereas the third uses experimental 
data. Finally, section 4 contains some conclusions. 



II. GENERAL THEORY 

For concreteness we consider a system with over- 
damped Langevin dynamics described by 



i = -V'{x,t) + ^/2j^£,{t) 



(2) 



where x denotes the degrees of freedom, V^ is a potential 
giving rise to a deterministic drift, and ^{t) is a stan- 
dard Gaussian white noise source obeying {^{t)) = and 
{^{t)£,{t')) = 5{t — t'). We denote derivatives with respect 
to a; by a prime and those with respect to t by a dot. 

During the time interval [0,ii] the potential changes 
between Vb(x) = V{x,0) and Vx{x) = V{x,ti) according 
to a fixed protocol. Using prepoint discretization the 
probability density functional of trajectories starting at 
t = at xo and ending a,t t = ti at xi is up to a constant 
given by 



p[a;(-)|a;o,a;i] ^ exp [ — P dt L{x{t),x{t),t) 



with the Lagrangian 

L{x,x,t)^^(^x + V'{x,t)y 



(4) 



The initial point, Xq, is sampled from the Gibbs measure 
corresponding to Vo{x) whereas the final point, xi, is free. 
For the work performed along a particular trajectory x{t) 
we have [Hi 



W[xi-)] = f'dtV{x{t),t) . 
Jo 

With the initial partition function 



Zo - dx e-^^«(^) 



(5) 



(6) 



(3) the probability distribution of the work is given by 



PiW) = — / dxoeM-l3Voixo)) / dx^ [ Vx{-)p[x{-)\xo,xi] SiW ^W[xi-)]) 



(7) 



2:(0)=2;o 



Using ^,^, and ^ we then find 



mn-j'^j...j 



dq 



x(ti)=xi 



Vx{-) e-0S[x(.)M 



(8) 



x{0)=X(, 



with the action 



S[x{-),q] = Fo(a;o) ^\j dt [\{x + V'f + iqV] - '-qW . 



(9) 



To apply the method of optimal fluctuations in the 
present context we evaluate the integrals in ^ for a pre- 
scribed value of W by the saddle-point approximation. 
Formally this corresponds to considering the weak noise 
limit /? — > oo. 



A. The optimal trajectory 

The determination of the optimal trajectory in ([5]) in- 
cludes the optimal choice of its initial and final point [3g| . 
Introducing the augmented Lagrangian 



L(x, i, t) = L(x, i, t) -I- i^V{x{t), t) (10) 



the corresponding Eulcr-Lagrange equation (ELE) takes 
the form 



d dL dL 



0. 



dt dx dx 
It is completed by the natural boundary condition 



(11) 



dL 

dx 



= 0, (12) 

at the end of the interval and the initial condition 

Uo'(xo) = , (13) 



dL 

dx 



incorporating the sampling of the starting point from the 
equilibrium distribution at t = 0. Solving pip - p3p and 



eliminating the Lagrange parameter q using ([5]) we gener- 
ically find for each value of W exactly one optimal tra- 
jectory x{t; W). The asymptotic estimate 



P{W) 



'l3S[x(-}M 



(14) 



for the distribution of work values becomes the more ac- 
curate the larger /3 is or, equivalently, the more W lies in 
thetailof F(M^). 



B. Neighbourhood of the optimal trajectory 

Although p^ gives a correct estimate of the asymp- 
totic behaviour of P{W) it is often desireable to im- 
prove its accuracy by incorporating the dominant pre- 
exponential factor. This factor has contributions from 
trajectories in the neighbourhood of the optimal one and 
also accounts for the Jacobian accompanying the transi- 
tion from p[a;(-)] to P(W). It is determined by including 
the quadratic fluctuations around the saddle-point into 
the calculation. This can be accomplished by adopting 
the Gelfand-Yaglom method i37, j38, 39] which yields an 
ordinary differential equation for the fluctuation deter- 
minant to the present problem. 

Two points are different from the standard case. First, 
the free endpoints of the optimal trajectory contribute 
to the Gaussian fluctuations and give rise to modified 
boundary conditions for the fluctuation determinant. 
Second, the constraint VF[x(-)] = W suppresses some 
fluctuations and gives rise to a correction factor to the 
free fluctuation determinant. Some details of the explicit 
calculation necessary to incorporate these two modifica- 
tions are given in the appendix. 

Using S = S[x{-),q], V = V{x{t),t) and similarly for 
derivatives of V the final result for the asymptotics of the 
work distribution is 



P{W) = 



-I3S 



= (l + 0(l//3)) (15) 



where Q{t) is the solution of the initial value problem 
= Q + 2V"Q + [(2 - iq)V" + (i - V')V"']Q 



Q{t = 0) = 1 
and R is given by 



Q{t = 0) = 



R^ I dt dt'V'{x(t),t) 
'o Jo 



S^S 



5x{t)5x{t') 



(16) 

V'(S(i'),t')- 
(17) 



III. EXAMPLES 



A. The shifted parabola 



As a first example we consider a Brownian particle 
dragged in a parabolic potential, i.e. 

V{x,t)^{x-tf/2. (18) 



This system has been been analyzed thoroughly both 
from the theoretical [40, |4l|, [43] as well as from the ex- 
perimental side 1431. The distribution P{W) is known to 
be Gaussian 140, SjJ 



P{W) 



27rcr^ 



exp 



{-, 



{W - a^/2) 



with 



'w 



2(ti 



1 



(19) 



(20) 



Since in this example the complete distribution P{W) is 
known exactly it merely serves as a test of our method. 
The ELE (fTTj) is for (fT8)) linear and can be solved an- 
alytically with the result 



x{t-W) = -(2i-t-e- 



M/^(2 - e-* - e*^*i 
2(fi+e-*i-l) 



This yields 



S[x{-),q\ = .,^ , ^_. 7T (21) 



4(ii 



1) 



which correctly reproduces the exponential factor in ([Tt 
The explicit form x{t] W) of the optimal trajectory for 
different values of W and ii characterizes the optimal 
combination of unlikely initial condition xq and rare re- 
alization of the noise ^(t) necessary to bring about large 
deviations in W . 

To determine the prefactor in (fT9|) we first observe that 
for (fT^ the differential equation p^ reduces to 

g + 2Q = o, Q(o) = i, g(o) = o (22) 

with the solution Q{t) = 1. Moreover 



(525, 



N 



Sx{t)dx{t') 



r.=-oS"it-t') + -d{t-t'). 



(23) 



Combining this expression with the boundary conditions 
2/(0) = j/(0) and y{ti) — —y{ti) for the fluctuations 
around the optimal path yields 



5^S 



N 



_5x{t)5x{t') 

With 1/' = -1 we then find 

i? = 2(ii-l- 



26{t-t')sinh{t-t') + e* 



e-^) 



'w 



(24) 



(25) 



Putting all together and using Zq — yj2'K/ (3 the prefac- 
tor of (fT9|) is also reproduced. In this simple example 
the asymptotic result hence already gives the complete 
distribution. 



B. The breathing parabola 

A more advanced example [3a] is given by the breath- 
ing parabola [4i.[45t. 



y(,,t)^M,2 



(26) 



for which the distribution of work is neither Gaussian nor 
completely accessible analytically. We will consider the 
case of a monotonously decreasing function k{t) implying 
W < and determine the asymptotic form of P{W) for 
VF -> -oo. The ELE ^ is given by 



+ {{l-iq)k~k^)x^O 



(27) 



whereas the boundary conditions (fT^ and (fT^ acquire 
the form 



i(0) = k{0)xQ 



and 



x{ti) = -k(ti)xi (28) 



respectively. These equations constitute a Sturm- 
Liouville eigenvalue problem which for the special choice 



m = 



1 



i+t 



(29) 



can be solved analytically. The result is 



J-W , 



V9UA 
hfi cos(/i ln(l + t)) + sin(^ ln(l + t))\ (30) 



where 



Km) = 2 (/^-— )sini^-cosi/-|-l+i^(Ai+— ) > 0. (31) 



and /i = \/iq — 9/4 is a solution of 



(V-3)sin- 



8u cos — = 
^ 2 



(32) 



with 1/ = 2/iln(l + ti). 

There are hence infinitely many discrete values 
qo,qi,... of q each associated with two trajectories 
x"^{t; W) and x^~{t\ W) related to each other by the in- 
version symmetry x — > —x of the problem. All x {t\ W) 
are local maxima oi p[x{-)]. However, it can be proved 
that p[a;''^(-)] > p[a;"^(-)] for all n > 0, i.e. the maxima 
at xP^{t\W) are the dominant ones. This is in accor- 
dance with intuition since large absolute values of W are 
realized by trajectories which are most of the time far 
from the minimum of the potential. On the other hand 
it is known from the general theory of Sturm-Liouville 
problems that the a;"^(i; W) have n zeros in the interval 
(0, ti). It is hence not surprising that the "ground state" 
solutions xP^{t\W) dominate the asymptotics of P(W^). 

Neglecting contributions from the sub-dominant max- 
ima we hence find from ([50]) . (pS)) . and ^ for the expo- 
nential term in the asymptotic work distribution 



P{W) 



-/3S[2"±(-),9o] 



JiKy.o)W 



(33) 




FIG. 1: Histogram of work values as obtained from 10^ simu- 
lations of the Langevin equation ((JJ with P = 2, ti — 100, and 
V{x,t) given by (|26|l and p9|) . The full line is the asymptotic 
form of the work distribution as derived in (|37p . The inset 
shows two estimates for the free energy difference AF together 
with their standard deviation as function of the sample size 
n. Circles are the standard estimate (|38|) . squares give the 
improved one (|39|) incorporating the asymptotic behaviour of 
P(W) with W* = —4/3. The dashed line is the exact result. 



where 



hi^^) 



8ff(A^) 



9 

(8^*2 - 6) cos ly - (2/1^ - 11/x -f — ) 

oil 



+ 8ii^ + 6 + iy{2ii^ + 5ii + —) 

8fi 



(34) 



Using (|3ip one can show that /i(/x) > 1 as is necessary 
for the existence of the Jarzynski average ^ . Note also 
that in the present case different values of W do not 
correspond to different values of q since the latter is fixed. 
As shown by (|30p different values of W are realized by 
different initial conditions of x^^{t; W). 

In the determination of the pre-exponential factor to 
the asymptotic we concentrate on its dependence on W. 
From (|16p we find using 



= Q + 2k{t)Q + {2 ^ iqo)k{t)Q (35) 

and therefore Q(ti) is independent of W. Likewise 



S^Sn 



5x{t)5x{t') 



l5"{t-t') + \{e 



{l-qo)k)S{t-t') 

(36) 
does not depend on W and hence neither does its inverse. 



On the other hand V' = kx°^ is proportional to \/—W 
as follows from ([50]) . This implies R ^ \f~W and we get 
the asymptotic result 



P{W) 



Jih(p.a)W 

4~w 



(37) 
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FIG. 2: Comparison of the simulation results (circles) for 
P{W) shown in Fig. [T] and the asymptotic behaviour (|37|) 
(full line) on a logarithmic scale. 



It is instructive to check this result against numerical 
simulations of the Langevin dynamics [3y|. Fig. [1] shows a 
histogram of work values obtained from such simulations 
together with the asymptotics ([57]) . Fig. [5] provides a 
logarithmic blowup of the small- ly region. To determine 
the prefactor in ([57]) a breakpoint W* is chosen and the 
area under the asymptotic form of P{W) for W < W* 
is set equal to the total weight $< of the histogram for 
W < W* (grey bars in FigH]). The value of W* has to 
be chosen such that on the one hand P{W) is already 
well approximated by its asymptotic form (|37p and on 
the other hand the region around W* is still sufficiently 
sampled by the histogram. As shown by Fig. [5] in the 
present case there is a whole window of admissible values 
of W* extending from roughly -1.5 down to around -5. 

The asymptotic form of P{W) can be utilized to im- 
prove the estimate ([T|) for the free energy difference AF. 
To show this we have subdivided the 10^ work values 
Wi obtained in the simulations into 10'' runs. Using 
n — 10^. ..10"^ values from each run we have then de- 
termined the standard Jarzynski estimate 



^^' -rM-n^^-""- 



(38) 



as well as an improved one 

W 



1 / r eC*-!)^^ 

AF™ = --ln(c / , 

P V J ^f-w 



n ^ — ' 



.-pw. 



n 

Wi>W' 



(39) 
using the asymptotic form of P{W) for W < W* . Here 
the constant c is determined from the normalization con- 
dition 



TV* 



$^ = c 



^hf3W 



V-w 



■.dW. 



(40) 



The inset in Fig[T] shows both estimates together with 
their standard deviation for different values of n as well 



as the exact result AF°''^'=* = -ln(l + ti)/{2(3). As is 
clearly seen both the bias and the standard deviation are 
significantly reduced when combining the histogram with 
the asymptotic form of P{W) as given by ([57)) . 



C. Driven Brownian particle near a wall 

We finally demonstrate the applicability of our method 
to the analysis of experimental data. In [3] a charged 
colloidal particle near a wall was subjected to a time- 
dependent anharmonic potential V{x^ t) generated by op- 
tical tweezers. Measuring the distance of the particle 
from the wall the distribution of work performed during 
one cycle of the potential modulation was determined 
(histogram in Fig. 4 in Q). Since Vq{x) = Vi{x) this case 
is characterized by AF = 0. 

As discussed in Q the dynamics of the particle may 
be approximately modeled by an overdamped Langevin 
equation. Due to the vicinity of the wall the friction 
coefficient and the noise intensity now depend on the 
state X. Moreover, in order to retain the Gibbs mea- 
sure as stationary distribution of the stochastic process 
an additional drift term has to be added [43|. Using Ito 
convention the resulting equation is [36| 



X = ~D{x)V'{x, t) + D'{x) + y/2D{x) ^{t) , (41) 
with potential 

V{x,t)=Ae-''^''-''^+B{t){x-a) (42) 

and state dependent diffusion coefficient [4^| 

Do 



D{x) 



R ■ 



(43) 



The values for Dq, the radius R of the particle, the pa- 
rameters A, K, a of V{x, t), and the protocol function B(t) 
are taken from the experiment [ig . Instead of ^ we now 
have 



L{x,x,t) ^ j^[x + D{x)V'{x,t) ^ D'{x)y 



(44) 



(45) 



)(x) 
The corresponding ELE 

D' 

^x - —x^ + {1 - iq)DV' 

+ {D' - DV'){DV" - D" + -D'V + — ) 

can no longer be solved analytically but its numerical 
solution does not pose any specific problems '36|. Solv- 
ing (|45p for a wide range of g-values and using the so- 
lution in ([5]) to establish the relation between q and W 
the extremal action S can be determined. The calcu- 
lation of the pre-exponential factor is now much more 
involved since the differential equation for Q{t) is more 
complicated and both Q{ti) and R will depend on W 
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FIG. 3: Histogram of 16200 work values obtained experimen- 
tally in [3| together with the asymptotic form (full line) de- 
rived from (|41|l - (|43| ). The inset shows the standard (circles) 
and the improved (squares) estimate for the free energy dif- 
ference. The exact result is AF = (dashed line) . 



in a non-trivial way. We leave this problem for further 
investigations and use for the present example only the 
asymptotic behaviour of P{W) resulting from p^ . It 
is shown in Fig. [3] together with the histogram of ex- 
perimental results. The inset gives again a comparison 
of the estimates for AF determined analogously to pS]) 
and dMI with W* = -3/2. 

While the asymptotic form of P{W) seems to be well 
captured the improvement in the free energy estimates is 
less distinctive than in Fig. [TJ The reason may be that 
P{W) decreases rather rapidly for small W which makes 
the matching between histogram and Asymptotics more 
difficult, in particular since the pre-exponential factor is 
not available. On the other hand, eq. (|4T|) is already 
an approximation to the experimental situation and the 
asymptotic form of P{W) derived from it may therefore 
differ from the true one. 



J 



IV. CONCLUSION 

We have shown that the method of optimal fluctua- 
tions allows to analytically characterize the asymptotic 
form of the work distribution in driven Langcvin systems. 
This information may be combined with histograms of 
work values as obtained in experiments or numerical sim- 
ulations to improve the accuracy of free energy estimates 
exploiting the Jarzynski equation. The method will work 
best in situations where an overlap region in W^- values ex- 
ists which is sufficiently sampled by the histogram and 
at the same time well described by the asymptotic be- 
haviour. 

Our method builds on a saddle-point calculation of 
a functional integral over stochastic trajectories con- 
strained to a specific value of the performed work W. Al- 
though similar techniques have been used in the context 
of non-equilibrium work and fluctuation theorems (see, 
e.g., [4^, Ho, [sil ) the application to constrained problems 
aiming at the asymptotic behaviour of the work distribu- 
tion is to our knowledge new. It will be interesting to 
generalize the method to higher-dimensional situations. 
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APPENDIX A 

In this appendix we give some details on the calcu- 
lation of the Gaussian fluctuations around the saddle- 
point in the integral ([5]). Using e = ti/N, tj = ej, 
Vj — V{x{tj),tj), and similar for the derivatives of 
V{x{t), t) the time-sliced version of this integral reads 



j=o 



P{W)^ lim -^(/-)f jdq l\\dx 
^ ' N^oo AttZq ^Ane' J J .J?, ' 

with the discretized action defined by 



PSN{{xj},q) 



(Al) 



N-l 



3=0 



\,w. 



(A2) 



Denoting the saddle-point values of Xj and q by an over- 
bar, using Sn = Sm{{xj}^ q) and expanding the exponent 



to second order in Xj — Xj and q — q we find 



P{W) = lim 



^ — Sn 



N-xx, Zq VdetM 



(l + 0(l//3)) 



(A3) 



where the symmetric matrix M is given by 



MkN+i = 2e 



'JV 



dxkdxi 



dxkdq 

'n 



d'^Sr. 



dqdq 



Aki for k,l = 0,...,N 



= ie^Vl for fc = 0, ...,7V 
= 0. 



Here Am is a tridiagonal fluctuation matrix of the usual 
form [33, [3^, [39|] . Its determinant can be obtained from 
a recursion relation which for e — > turns into a dif- 
ferential equation. Analogous to the standard Gelfand- 
Yaglom procedure we find det A = Q{ti) where Q{t) is 
the solution of the initial value problem (fT6|) . 

In order to reduce the calculation of det M to that 
of det A we multiply the first iV + 1 rows of M by 
—ie^{V')'^A~^ and add this to the last row. The re- 
sulting matrix has then in the last row all zeros except 



for the last entry which reads 



(A4) 



k,l 



Consequently det M = Rn det A. This result is in fact 
quite intuitive. Assume for simplicity that the constraint 
MK[a:(-)] = W is orthogonal to one eigenvector e„ of A 

with eigenvalue A„. Then {Vl} which is the gradient of 
the constraint is parallel to e„ and Rn is proportional to 
1/A„. It hence cancels exactly that eigenvalue of the un- 
constrained fluctuation matrix A describing fluctuations 
perpendicular to the constraint which are forbidden. 
Using 



AT- 



lim Rm — e R 



(A5) 



with R given by (|17p we finally end up with (fT5 
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